home *** CD-ROM | disk | FTP | other *** search
/ Libris Britannia 4 / science library(b).zip / science library(b) / MATH / FSULTRA1.ZIP / ULTRAMAC.C < prev    next >
Text File  |  1992-06-18  |  7KB  |  208 lines
  1. /* 
  2. FSU - ULTRA    The greatest random number generator that ever was
  3.         or ever will be.  Way beyond Super-Duper.
  4.         (Just kidding, but we think its a good one.)
  5.  
  6. Authors:    Arif Zaman (arif@stat.fsu.edu) and
  7.         George Marsaglia (geo@stat.fsu.edu).
  8.  
  9. Date:        27 May 1992
  10.  
  11. Version:    1.05
  12.  
  13. Copyright:    To obtain permission to incorporate this program into
  14.         any commercial product, please contact the authors at
  15.         the e-mail address given above or at
  16.  
  17.         Department of Statistics and
  18.         Supercomputer Computations Research Institute
  19.         Florida State University
  20.         Tallahassee, FL 32306.
  21.  
  22. See Also:    README        for a brief description
  23.         ULTRA.DOC    for a detailed description
  24.  
  25. -----------------------------------------------------------------------
  26. */ 
  27. /*
  28.    File: ULTRA.C
  29.  
  30.    This is the ULTRA random number generator written entirely in C.
  31.  
  32.    This may serve as a model for an assembler version of this routine.
  33.    The programmer should avoid simply duplicating and instead use the
  34.    usual assembler features to increase the speed of this routine.
  35.  
  36.    Especially the subroutine SWB should be replaced by the one
  37.    machine instruction (usually called subtract-with-borrow) that
  38.    is available in almost every hardware.
  39.  
  40.    For people not familiar with 8086 assembler, it may help to
  41.    consult this when reading the assembler code. This program should
  42.    be a dropin replacement for the assembler versions, but is about
  43.    half as fast.
  44. */
  45.  
  46. #define N  37           /* size of table        */
  47. #define N2 24           /* The shorter lag      */
  48.  
  49. long swb32[N],        /* These arrays contian the seed mixed with    */
  50.      swb16[N],         /* a congruential. They are used 32, 16, and 8  */
  51.      swb8[N],         /* bits at a time respectively.            */
  52.  
  53.      swbseed[N] = { 0x0D45D69A,    /* The seed for the SWB generator    */
  54.   0x9DB73B1A, 0xA84604E8, 0x7C5F0CA5, 0xBC0196CE, 0xFF4CB42E, 0x7ACA6BE3,
  55.   0xA9ED2A5A, 0x6405A8F7, 0xAC00D4F8, 0xBDD1FC77, 0x064F9DC5, 0xF10AB737,
  56.   0x781293D3, 0x8410C2D2, 0x1C6587DB, 0x7D8F8F0F, 0xF3DCC4EA, 0xB965F99F,
  57.   0x9A0094D1, 0x65976D1C, 0x09173DC1, 0x8E38B992, 0x84701D44, 0x14F0E1B9,
  58.   0xE8A1EC5F, 0x1E925A12, 0xE77A0B5B, 0xCDB5926E, 0xD16260C8, 0xC917E806,
  59.   0x519076AA, 0x7BF6C21C, 0x808C0C90, 0x3E93C7B8, 0x707D6EA0, 0xF1DB698D};
  60. char    flags=0;        /* The carry flag for the SWB generator    */
  61. long unsigned congx = 0x1C3a128F;    /* Seed for congruential generator    */
  62.  
  63. long    *swb32p;    short    swb32n=0;    /* A counter and a    */
  64. short    *swb16p;    short    swb16n=0;    /* pointer is kept for    */
  65. char    *swb8p;        short    swb8n=0;    /* each array.        */
  66.  
  67. char     swb1[32];                /* The swb1 is slightly    */
  68. char    *swb1p;        short    swb1n=0;    /* different.        */
  69.  
  70. /*************************************************************************
  71.  * For each of swb32, swb16, swb8 and swb1, there is an array, a counter *
  72.  * and a pointer. The counters count down, the pointers go up. When the  *
  73.  * counter reaches zero, a fill routine is called to refill the array.   *
  74.  ************************************************************************/
  75.  
  76. /* rinit initializes the constants and fills the seed array one bit at
  77.    a time by taking the leading bit of the xor of a shift register
  78.    and a congruential sequence. The same congruential generator continues
  79.    to be used as a mixing generator for the Subtract-with-borrow generator
  80.    to produce the `ultra' random numbers
  81.  
  82.    Since this is called just once, speed doesn't matter much and it might
  83.    be fine to leave this subroutine coded just as it is.
  84.  
  85.    PS:    there are quick and easy ways to fill this, but since random number
  86.     generators are really "randomness amplifiers", it is important to
  87.     start off on the right foot. This is why we take such care here.
  88. */
  89.  
  90. void rinit(long unsigned congy,long unsigned shrgx)
  91. { short i,j;
  92.   unsigned long  tidbits;
  93.   congx=congy*2+1;
  94.   for (i=0;i<N;i++)
  95.   { for (j=32;j>0;j--)
  96.     { congx = congx * 69069;
  97.       shrgx = shrgx ^ (shrgx >> 15);
  98.       shrgx = shrgx ^ (shrgx << 17);
  99.       tidbits = (tidbits>>1) | (0x80000000 & (congx^shrgx));
  100.     }
  101.     swbseed[i] = tidbits;
  102.   }
  103.   swb32n = swb16n = swb8n = swb1n = 0;
  104.   flags = 0;
  105. }
  106.  
  107. /* SWB is the subtract-with-borrow operation which should be one line
  108.    in assembler code. This should be done by using the hardware s-w-b
  109.    operation in the SWBfill routine.
  110.  
  111.    What has been done here is to look at the msb of x, y and z=x-y-c.
  112.    Using these three bits, one can determine if a borrow bit is needed
  113.    or not according to the following table:
  114.  
  115.     msbz=0  msby=0  msby=1          msbz=1  msby=0  msby=1
  116.  
  117.     msbx=0  0       1               msbx=0  1       1
  118.     msbx=1  0       0               msbx=1  0       1
  119.  
  120.    PS: note that the definition is very carefully written because the
  121.    calls to SWB have y and z as the same memory location, so y must
  122.    be tested before z is assigned a value.
  123. */
  124. #define SWB(c,x,y,z) \
  125. c = (y<0) ? (((z=x-y-c) < 0) || (x>=0)) : (((z=x-y-c) < 0) && (x>=0));
  126.  
  127. void SWBfill(long *x)
  128. { short i;
  129. /*
  130.    The following two lines are the heart of the system and should be
  131.    written is assembler to be as fast as possible. It may even make sense
  132.    to unravel the loop and simply wirte 37 consecutive SWB operations!
  133. */
  134.   for (i=0;  i<N2; i++) SWB(flags,swbseed[i+N-N2],swbseed[i],swbseed[i]);
  135.   for (i=N2; i<N;  i++) SWB(flags,swbseed[i  -N2],swbseed[i],swbseed[i]);
  136. #ifndef fast
  137.   for (i=0;  i<N;  i++) *(x++) = swbseed[i] ^ (congx = congx * 69069);
  138. #else
  139.   for (i=0;  i<N;  i++) *(x++) = swbseed[i];
  140. #endif
  141. }
  142.  
  143. long swb32fill(void)
  144. { swb32p = &swb32[0];
  145.   SWBfill(&swb32[0]);
  146.   swb32n = N-1;
  147.   return *(swb32p++);
  148. }
  149.  
  150. short swb16fill(void)
  151. { swb16p = (short *) &swb16[0];
  152.   SWBfill((long *) swb16p);
  153.   swb16n = 2*N-1;
  154.   return *(swb16p++);
  155. }
  156.  
  157. char swb8fill(void)
  158. { swb8p = (char *) &swb8[0];
  159.   SWBfill((long *) swb8p);
  160.   swb8n = 4*N-1;
  161.   return *(swb8p++);
  162. }
  163.  
  164. #define i32bit()  ((swb32n--) ? *(swb32p++) : swb32fill())
  165. #define i31bit() (((swb32n--) ? *(swb32p++) : swb32fill()) & 0x7FFFFFFF)
  166. #define i16bit()  ((swb16n--) ? *(swb16p++) : swb16fill())
  167. #define i15bit() (((swb16n--) ? *(swb16p++) : swb16fill()) & 0x7FFF)
  168. #define i8bit()   ((swb8n--)  ? *(swb8p++)  : swb8fill() )
  169. #define i7bit()  (((swb8n--)  ? *(swb8p++)  : swb8fill() ) & 0x7F)
  170.  
  171. char swb1fill(void)
  172. { unsigned long i,j;
  173.   swb1p = &swb1[0];
  174.   swb1n = 31;
  175.   i = i32bit();
  176.   for (j=0;j<32;j++) { swb1[j] = i&1; i=i>>1; }
  177.   return *(swb1p++);
  178. }
  179.  
  180. #define i1bit()   ((swb1n--)  ? *(swb1p++)  : swb1fill() )
  181.  
  182. #define two2neg31  ( (2.0/0x10000) / 0x10000 )
  183. #define two2neg32  ( (1.0/0x10000) / 0x10000 )
  184.  
  185. float vni(void)
  186. { long temp;
  187.   temp = i32bit();
  188.   if (temp & 0xFF000000) { return temp * two2neg31; }
  189.   return (temp + i32bit() * two2neg32) * two2neg31;
  190. }
  191.  
  192. float uni(void)
  193. { long temp;
  194.   temp = i31bit();
  195.   if (temp & 0xFF000000) { return temp * two2neg31; }
  196.   return (temp + i32bit() * two2neg32) * two2neg31;
  197. }
  198.  
  199. double dvni(void)
  200. { unsigned long temp;
  201.   temp = i32bit();
  202.   return ( (long) i32bit() + temp * two2neg32) * two2neg31; }
  203.  
  204. double duni(void)
  205. { unsigned long temp;
  206.   temp = i32bit();
  207.   return ( i31bit() + temp * two2neg32) * two2neg31; }
  208.